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The viscosity and self-diffusion constant of particle-based mesoscale hydrodynamic methods, 
multi-particle collision dynamics (MFC) and dissipative particle dynamics (DPD), are investigated, 
both with and without angular-momentum conservation. Analytical results are derived for fluids 
with an ideal-gas equation of state and a finite-time-step dynamics, and compared with simulation 
data. In particular, the viscosity is derived in a general form for all variants of the MFC method. 
In general, very good agreement between theory and simulations is obtained. 
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I. INTRODUCTION 

Soft matter systems such as polymer solutions, col- 
loidal suspensions, membranes, and microemulsions ex- 
hibit many interesting dynamical behaviors, where hy- 
drodynamic flow plays an important role, as do thermal 
fluctuations. The characteristic time and length scales 
of soft-matter systems are in the range from nanosec- 
onds to seconds and from nano- to micrometers, respec- 
tively, and are thus typically much larger than the atom- 
istic scales. Mesoscale simulation techniques are there- 
fore necessary to simulate these systems for sufficiently 
large system sizes with reasonable computational effort. 
Several mesoscale techniques for the simulation of the 
flow of complex fluids accompanied by thermal fluctua- 
tions have been developed in the last decades, such as 
direct simulation Monte Carlo (DSMC) P, the Lat- 
tice Boltzmann method 3, 4|, dissipative particle dy- 
namics (DPD) a i 0, i 1 M El, M M M M 
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DPD, and MPC are off-lattice hydrodynamics meth- 
ods and share many properties. DPD and MPC have 
been applied to various soft-matter systems such as col- 
loids [iHl, M, [m, [Sl^olymers Jl ijEl [13, [H, Hi, IH , 

and surfactants 

The key features to distinguish DPD and MPC are the 
application of a Langevin thermostat to the relative ve- 
locities of particle pairs or multi-particle collisions, and 
whether or not to employ collision cells. To understand 
and elucidate the relation between DPD and MPC, two 
intermediate methods have been proposed in Ref. [20j . 
which are DPD with a multibody thermostat (DPD-MT) 
and MPC-Langevin dynamics (MPC-LD). The standard 
MPC algorithm does not conserve angular momentum. 
However, an angular-momentum-conserving version of 
MPC has also been proposed in Ref. [2(|. We denote the 
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versions of a simulation method with or without angular- 
momentum conservation by an extension '-fa' or '—a', 
respectively. The importance of angular-momentum con- 
servation in MPC fluids has been studied in Ref. ^2[. In 
the absence of angular-momentum conservation, an ad- 
ditional torque appears which depends linearly on the 
vorticity, whereas the velocity field is unaffected. There- 
fore, it is essential to employ '-|-a' techniques to simulate 
systems such as rotating colloids and binary fluids with 
different viscosities. 



In this paper, we investigate the viscosity rj and self- 
diffusion constant D of MPC and DPD methods. The 
transport coefficients of '—a' versions of MPC were pre- 
viously derived analytically, and show good agreements 
with numerical results [M, HE 113, [11,11^. We derive 
here analytically the viscosity and diffusion constant of 
all '+a' versions of MPC. 



The transport coefficients of original version of DPD 
were derived analytically for systems with an ideal-gas 
equation of state in the small-time-step limit [T^ and 
with finite time step (2l| . and phenomenologically for 
soft-repulsive interactions [2l|. Here, we investigate 
the transport coefficients of DPD— a and DPD-MT for 
the ideal-gas equation of state with finite time step. 
The viscosity and diffusion constant are also determined 
from simulations of simple shear fiow with Lees-Edwards 
boundary conditions and of the mean square displace- 
ment of a particle, respectively. 



The outline of this paper is as follows. In Sec. [TTl 
we describe several versions of MPC, both with and 
without angular momentum conservation, and calculate 
their transport coefficients analytically and numerically. 
Transport coefficients of several version of DPD are cal- 
culated in Sec, mil In Sec. lIVi we discuss the upper hmits 
of the local shear rate for which thermostats in MPC and 
DPD are capable to provide local-equilibrium condition. 
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II. MULTI-PARTICLE COLLISION DYNAMICS 
(MPC) 

A. Simulation Method 

1. MPC without angular-momentum conservation 

MPC is a modification of DSMC to include multi- 
particle collisions, in order to make the algorithm more 
efficient in its application [l^]- A fluid is described by 
point-like particles of mass m. The MPC algorithm con- 
sists of alternating streaming and collision steps. In the 
streaming step, the particles move ballistically, 



(1) 



where Ai is the time interval between collisions. In the 
collision step, the particles are sorted into cubic cells of 
lattice constant /c- The collision procedure is different for 
each version of MPC. For MPC— a, it is generally given 
by 



(2) 



where is the velocity of the center of mass of all par- 



ticles in the box, and 



\i — v^. The collision op- 



erator r2[vi_c] stochastically changes the relative velocity 
Vi^c, with X^ieccii ^[^i-c] = to keep the translational 
momentum constant. This stochastic process is indepen- 
dent for each cell and each time step, and the collision 
operator r2[vi^c] depends on whether a particle is inside a 
cell, but not on its position within the cell. To guaran- 
tee isotropy, the operator must be symmetric on average, 
with {va^\y]p) = (1 — A){va^)5af3^ where the subscripts 
a, l3 G {x, y, z\ indicate the spatial components. The 
constants A and B = \ — {Vl[v\a^[v\p) / {vaVp) are char- 
acteristic quantities of each version (see Table H]), which 
play an essential role in determining the transport coef- 
ficients. The operator J7[vj^c] conserves the total kinetic 
energy in each cell (local micro-canonical ensemble) or is 
coupled to a thermostat (local canonical ensemble) . The 
collision cells are randomly shifted before each collision 
step to ensure Galilean invariance [23j . 

The operator r2[v] of the original version of MPC is 
the rotation operator. It is represented by a matrix 
J^r(v) which rotates velocities by an angle 9. The ro- 
tation axis is chosen randomly for each cell, which re- 
quires one integer or two real random numbers in two- 
(2D) or three-dimensional (3D) space, respectively. In 
2D, the axis is the ±z direction (out of plane), i.e the 
rotation is clockwise or anticlockwise with the angle 6 
(see Fig. [1]). This original version of MPC is typically 
denoted MPC or stochastic rotation dynamics (SRD). 
We denote it MPC-SR— a in this paper, in order to dis- 
tinguish this particular version clearly from the whole 
family of MPC techniques. In MPC-SR— a, the energy 
in each cell is conserved. The temperature can be con- 
trolled by an additional rescaling of the relative velocities 



spatial dimension, N is the total number of particles, and 
iVccii is the number of cells occupied by particles. This 
corresponds to a velocity-sca ling version of the profile- 
unbiased thermostat (PUT) [44!|, where cells are intro- 
duced to thermostat local velocities relative to the center- 
of-mass velocity of each cell. The number d{N — A^coii) 
of the degrees of freedom should be sufficiently large for 
the central-limit theorem to apply. This usually implies 
that the number of cells included in the calculation of the 
rescaling factor is large. When the velocity rescaling is 
performed on the level of single collision cells, the Monte 
Carlo scheme proposed in Ref. [S^l should be employed. 

In the random angle version of MPC (denoted MPC- 
RA— a) [2^, the same matrix J1r(v) is employed, but the 
rotational angle 9 is also selected stochastically varied in 
the interval < 9 < 9o- In MPC-RA-a, one or three 
real random numbers are required for each cell in 2D or 
3D, respectively. 

In the Andersen-thermostat [isl . \^ version of MPC, 
denoted MPC- AT [13, HI], the operator completely re- 
news the relative velocities in the cell, J7[v] — v™" — 
'^5^"/A^c, where Nc is the number of particles in a 
cell. A velocity vj'^" is chosen from a Maxwell-Boltzmann 
distribution. Thus, in MPC- AT— a, the velocities of par- 
ticles are updated by 



J G cell 



J 



(3) 



Instead of the energy, the temperature is constant in 
MPC- AT. 

In the Langevin version of MPC (MPC-LD-a) 0, 
the Langevin thermostat is applied to the relative veloc- 
ities in a collision cell. The particle motion is governed 
by 



dv, dU 

m— = 7Vi,c + 

at OTi 



J G cell 



Nr. /■ 



(4) 



In order to satisfy the fluctuation-dissipation theorem, 
the Gaussian white noise ^j(i) has to have the aver- 
age (^i,a(t)) = and the variance {£,i,ait)£,j,p{t')) = 
2kBT5ij5ap5(t — t'), where a,/3 G {x,y,z} and fcer is 
the thermal energy. We consider in this paper only fluids 
with an ideal-gas equation state, i.e. t/ = in Eq. Q. 
The finite time-step version of MPC-LD— a is given by 
the leapfrog algorithm, 

r»(^n+l/2) = r»(^ri-l/2) + Vj,„At, (5) 

Vj(i„+i) = v^+aidVj,c(in) + 6id|Cj,„- E %^}' 



l-7Ai/2m 
with aid = :;— , Old 



l + -fAt/2r 



bid = 



l + -fAt/2r 



(6) 



Nccii)kBT/mJ2i^-^ 



j^c^, where d is the 



where {^i,n,a) = and {^r,n,a^j,n' ,p) = '^kBTdijSaflSnn' 

Thus, the collision operator is r2[vi^c] — aidV^.c 4 
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TABLE I: Correlation factors A = \~ {va_^\v\a) I (va^) and 
B = \ — (^\\r\a^\w\(j) I {vaVp) ofvarious MPC methods, where 
a, /3 £ \x, y, z} and a ^ p. 



A 



B 
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2/1 sin So \ 




1 sin 2Bq 
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r) sin ^0 sin 2Bq \ 
Bo 200 ' 


(d = 
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MPC-AT 


1 




1 






MPC-LD 


^Ai/ m 




2^At/m 






l + 7At/2m 




(l + 7At/2m)^ 








-a 




old 



+a-vs 



DR 



+a-vs 
+a 



bld{i^ n-J2^3 n/Nc}- MPC-LD with jAt/2m = 1 coin- 
cides with MPC-AT. In MPC-AT and MPC-LD, the cor- 
relations have a simple relation, (\ — B) — (\ — Af' . How- 
ever, MPC-SR and MPC-RA have additional correlations 
between x and y components, i.e. (1 — i3) ^ (1 — Af' as 
shown in Table [J 



FIG. 1: (Color online) Schematic representation of the col- 
lision operation for MPC-SR±o and MPC-DR in 2D in the 
co-moving reference frame (with J^Vi = 0) at A'c = 3 and 
Q — n/2. Circles represent the positions of particles (o) and 
the center of mass (•). 'old' indicates the velocities before the 
collision, '±a' and 'DR' represent the velocities after the colli- 
sion for MPC-SR±a and MPC-DR, respectively, and '-|-a-vs' 
indicates the velocities after the '-|-a' collision without veloc- 
ity rescaling. 



2. MPC with angular-momentum conservation 

Collisions described by Eq. ([2]) conserve translational 
momentum, but do not conserve angular momentum. 
However, angular-momentum conservation can be im- 
posed by an additional constraint. This modification 
is straightforward for the MPC versions with an intrin- 
sic thermostat (such as MPC-AT and MPC-LD). In this 
case, the collision is given by 



mil 



J G cell 



(7) 



{rj,c X (vj,c - rJ[vj,c])} 



where II is the moment-of-inertia tensor of the particles 
in the cell. The relative position is r^^c = — where 
is the center of mass of the particles in the cell. The 
angular momentum of the cell after the collision, Ilwc = 
m J2 ^j.c X c , is the same as before the collision. The 
subtraction of either position or velocity of the center of 
mass can be omitted in the last term of Eq. ([7]), since 

E I-i,c X Vj- c - E X Vj- c = E I-i,c X Vj. 

For MPC-AT+a or MPC-LD+a, the terms 



fAT+a = mn-i J2 {""J^- ^ (vj - vf")} X r,,e, (8) 
fLD+a = mU'^ |rj,c X {jVj - ^/7lj(t)}| x r,,c 

(9) 



are added to Eqs. ^ and (|4]), respectively [20j . 

When Eq. ^ is apphed to the operator of MPC-SR 
or MPC-RA, the kinetic energy is not conserved. Thus, 
the collision process has to be modified by combining it 
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FIG. 2: (Color online) Radial distribution function g(r) of 
MPC-SR-fa (with n = 1, At* = 1 and n = 10, At* = 0.1) 
and MPC-AT -l-a (with n = 1, At* = 0.1) in two-dimensional 
space. The inset shows the n dependence of gir) of MPC- 
SR-|-a at At* = 0.1. Error bars are shown at several data 
points. 



with velocity rescaling to conserve the energy, 

' mUr^ ^ (r^^c X Vj,c) X r^^c (10) 



J G cell 



-I- 0|rj[vi,c] - mil ^ ^ (rj,c X Jl[vj,c]) X ri,c|, 

J G cell 

where <\> = {E,econ«)n/{E,eccn(u?)'}- Here, 
the relative velocities before and after the collision, 
u°'^ and u^, respectively, are given by = v^^c — 
toII^^ Ejecou(''i:C X Vj^c) X ri_c, where the total trans- 
lational and angular velocities of the cell are subtracted. 
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This collision is shown schematically in Fig. [TJ Under 
the molecular-chaos assumption, this yields the ideal-gas 
equation of state. However, the molecular- chaos assump- 
tion is not perfectly valid. Thus, the radial distribution 
function g{r) of MPC-SR+a exhibits deviations from the 
uniform distribution of the ideal gas, in particular for 
small n and small At (see Fig. [5]) . If the velocity rescaling 
for the energy conservation is done not for each cell but 
for the sum of many cells, this deviation becomes larger. 
A similar deviation is seen in DPD simulations Q with 
the modified velocity- Verlet algorithm MFC- AT -Ha 
and MPC-LD-ha and all '-a' versions of MFC give the 
correct uniform g{r) — see, e.g., the data of MFC-AT-|-a 
in Fig.d Thus, MFC-SR-ha should not be used for smaU 
n or small At. We recommend to check g(r) for any new 
MFC operator. 

An alternative modification of MFC-SR for two- 
dimensional fluids to conserve angular momentum has 
been proposed recently by Ryder (431 (see also Ref. [ist). 
We denote this algorithm MFC-DR (deterministic rota- 
tion). In MFC-DR, a rotational angle is chosen determin- 
istically to keep the total angular momentum of particles 
in a collision cell constant by the requirement W{1 — 
cos(6')} -I- Q sin(6') = 0, where W = X^jgcgH ^j.c x '^j-c and 
Q = EjGceii ""j.c ■ '^i.c- This implies 

W^^Q^ . 2WQ 

The velocities after a collision in MFC-DR are different 
from those in MFC-SR-|-a, since the '-|-a' procedure (from 
'—a' to '-fa-vs' in Fig.[T]) does not change the radial ve- 
locities. MFC-DR gives the correct uniform g{r) and is 
less time-consuming than other '-|-a' versions of MFC. 
We also checked that MFC-DR yields the correct con- 
stant angular velocities for phase-separated binary fluids 
with different viscosities in a circular Couette flow, as de- 
scribed in Sec. IV. C of Ref. ^4^]. However, this algorithm 
cannot be generalized to three-dimensional systems. 

B. Transport Coefficients 

1. Stress tensor 

Angular-momentum conservation implies that the 
stress tensor aaf3 for an isotropic Newtonian fluid is sym- 
metric, i.e. Uafi — (7/3a [48]. In contrast, MFC— a fluids 
have an asymmetric stress tensor 

(Jaf! = X{V-^)S^p (12) 

because of the lack of angular-momentum conservation 
[27L [29I . l4l |. where a,(3 & {x,y, z} and A is the second 
viscosity coefficient, f] and 77 are the symmetric and anti- 
symmetric components of the viscosity, respectively. The 



last term in Eq. ([12]) implies that the stress depends lin- 
early on the vorticity V x v, and does not conserve an- 
gular momentum. Thus, this term vanishes {i.e. = 0) 
in angular- momentum conserving systems. 

The evolution of the velocity field v(r) is determined 

by 

= -VP + (A + 7? - 7?) V(V • v) + (^ + 77) V^v, (13) 

where D/Dt is Lagrange's derivative and P is the pres- 
sure field. When a fiuid is incompressible, Eq. (fT3|) is the 
normal Navier-Stokes equation with viscosity t] ~ fj + fj. 
This is consistent with the usual definition of the shear 
viscosity rj = a^y/j in simple shear flow with velocity 
field V — 77/63;, where is the unit vector along the 
X direction. Since both the equations of continuity and 
of velocity evolution are of the same forms in systems 
with and without angular-momentum conservation, the 
absence of angular-momentum conservation does not af- 
fect the velocity field of a fiuid when the boundary con- 
ditions are given by velocities. However, it generates an 
additional torque, as described in detail in Ref. [HI. In 
this paper, we discuss the stress tensor of various MFC 
and DFD methods. 



2. MPC wittiout angular-momentum conservation 

The shear viscosity is calculated from Uxy/i = V = 
77 + 77 in simple shear fiow with v — 'jyex . The viscosity 
of MFC fiuids consists of two contributions, 77 = 77kin + 
77coi, where the kinetic viscosity 77kin and the collisional 
viscosity rjcoi result from the momentum transfer due to 
particle displacements and collisions, respectively. The 
derivation of the viscosity for MFC-SR— a described in 
Refs. [2^, [13, [H, [2^ can be employed directly for other 
'—a' versions of MFC, since the differences appear only 
in the factors A and B listed in Table [J 

The kinetic stress cr^JJ^ = 77i;in7 is the momentum fiux 
due to particles crossing a xz plane at ?/ = t/q = 0. The 
stress due to streaming in the time interval [t,t + At] is 
written as 

kin ^ f \ ^ 

''^y ^ jAti ^ 

ydt)<0,v,,y>-^ 

where 5* is the surface area of the considered plane. The 
average over equivalent xz planes yields 

f^x'y = -mn{VxVy)t+^t/2 = ~y Vi^xVt,y, (15) 

i 

where n ~ (Nc) is the average number of particles per 
cell, and V is the volume of the considered region V, 
with r-i e V; here the middle position ri{t -\- At/2) = 



5 



ri(t)+ViAt/2 during streaming is employed to determine 
whether the i-th particle is inside the region V. The ex- 
pression (|15p is symmetric in x and y. The symmetry 
of the kinetic part of the stress tensor, i.e. a^™ = cr^J)\ 
implies Tykin — for all versions of MPC and DPD. Nu- 
merically, a^y^ and jykin can be calculated from Eq. p4|) 
or The velocity distribution is shifted by particle 

streaming so that 

{vxVy)t.t+At = J dw VxVyPy{v^jvyAtex/2) 

= {VxVy)t+At/2 ± {vy'^)jAt/2, (16) 

where Pvi'v) is the velocity probability distribution. The 
velocity distribution is modified by the MPC collisions 
so that {vx^'^Vy'''''") = (1 - c,^){vxVy), where the factor 
Cm is determined later. The self-consistency condition of 
a stationary shear flow is {vxVy)t = {vxVy)t+At = (1 ^ 
Cm)i(vT. t ~7 At (ji,, ^ ) ) . The kinetic viscosity ?7ki„ is then 
given by [26| 

_ nk^TM / I 1^ 

Eq. (dZD holds for aU '±a' versions of MPC and DPD. 

The velocity correlations for MPC— a are calculated by 
using Eq. ©, 

lA + -fl--)(l-^) 



(17) 



/..new ncw\ 
\^i,x i,y I 



+- 



1 n2 

2A-B 



1 



{^-B(l-—)]{v,^xV^,y), (18) 



where molecular chaos is assumed, 

^^i,x'^j.y 

tion factor for a cell occupied by A^c particles is c{Nc) — 
B(l — 1/Nc). An MPC fluid is thermodynamically an 
ideal gas, so that the cell occupation number Nc fluctu- 
ates with the Poisson distribution, P{Nc) = e~"n^'' /Ncl 
with n — (Nc). Thus, the average correlation is give 
by Cm = Efeli c{k)P{k)k/n ^ B{n - I + e-")/n. The 
kinetic viscosity of MPC— a is then given by 

nkBTAtr n/B 
I n 



^i^x^i^yl 

Vj^xVj^y) and (vi^xVj.y) — for i ^ j. Thus the correla- 



(19) 



The collisional stress cr^°' — ?7coi7 is the momentum 
flux due to MPC collisions in cells crossing a plane at 
y = yo- It is given by ^] 

_col /,,ni 



Ir 



(20) 



When Eq. ((20|) is averaged over the planes crossing the 
cell, j/cc — ^c/2 < yo < Vcc + the stress reads 

E (t^ + I) «r (21) 



_col 



'xy 



Ic. 



'At 



where j/i^cc = Vi — Vcc and j/cc is the j/ component of the 
center-of-cell position Tcc- Numerically, g™ ^ and 77coi can 
be calculated from either Eq. ([20| or ([21]) . The mean 
velocity difference is (w,'^™ - Vi^x) = -(1 - ■^)A-jyi^cc, 
because (v^) — Vi/Nc, where j/j is averaged over —lc/2 < 
yj < lc/2 for j 7^ i at y^c = 0. Then the collisional 
viscosity rycoi of MPC— a is given by 



??col 



{ 5] (iVe - l)P(7Ve) 



Am 

Am{n — 1 + e^") 



12Z, 



d-2 



At 



(22) 



The vorticity viscosity 77C017 = 



col _ col 
xy yx 



is proportional 



to the angular-momentum transfer with respect to the 
origin {xcc - ^c/2,ycc - ^c/2), see Eq. Thus, the 

vorticity viscosity of MPC-l-a vanishes, 77C0I = 0, be- 
cause of angular-momentum conservation. For MPC— a, 
the molecular-chaos assumption gives a™' = 0, because 
{Vy™{x)) = {vyix)) = 0. Thus, the viscosities are 
V = Vcoi — ?7coi/2 [23,14^. This viscosity relation holds 
for all '—a' versions of MPC and DPD described in this 
paper. 

As an extension of this approach, the angular- 
momentum constraint can be applied only partially, by 
employing alternatively the MPC-coUision algorithms 
which conserve [given by Eq. ([7|)] and do not conserve 
[determined by the difference of the right-hand sides of 
Eqs. (m and ([7])] angular momentum. In this way, the 
viscosity ratio fi/rj can be varied continuously between 
and approximately 1. 

Next, we derive the self-diffusion constant D of 
MPC— a. Under the molecular-chaos assumption, 
the velocity correlation function decays exponentially, 



Vi^x{kAt)vi^xiO)) 



'-Ji,x 

{v^:•v^,x 

hym 



= (1 " Sm) kBT/m with 1 — Sm = 
The diffusion constant is thus given 



_ ksTAt / 1 
m Vs„, 



k=l 



:,,(fcAt)v,,,(0))} 



(23) 



IS Sn 



In MPC— a, the correlation factor 

s(fc)P(fc)fc/7i = yl(n - 1 + e-")/n with 
s{Nc) = A{1 - 1/iVc); this implies 



D 



k^TAt ( n/A 



\n- 1 



(24) 



However, the velocity auto-correlation func- 
tion {vx{kAt)vx{Q)) for small mean free path 
l\ = At-^/fceT/mo has a long-time tail due to hy- 
drodynamic backflow [1^, [3l|, |3^. This leads to an 
additional hydrodynamic contribution to the diffusion 
constant _D, which thereby becomes larger than predicted 
by Eq. 
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3. MFC with angular-momentum conservation 

To derive expressions for the self-diffusion constant and 
viscosity of MPC+a, we employ Eqs. ([T7D, and 
(I23p . which remain valid with angular-momentum con- 
servation. However, the correlation factors Sm and Cm 
of MPC+a are different from those of MPC-a. First, 
we consider the limit of large n, where Sm = and 
Cm — c(n), and derive the corrections for small n sub- 
sequently. The velocity correlation is calculated from 



X V, 



with the molecular-chaos assumption. The positions of 
particles are averaged over the cell, so that r^^c^ = 
(1 - l/iVc)/c^d/12 and H = (iV^ - l)ml^^I/6 where I 
is the identity matrix. Angular-momentum conservation 
implies additional correlations, which result in 



A I 



1 



2Nc 



(25) 




where Xi^cc is the x component of unit vector fi^cc — 
ri,cc/7'i,cc and (if cc) — V*^- The diffusion constant of 
MPC+a for large n is thus found to be 



D 



kBTAt f 



n/A 



\n- l)/2 



(26) 



For the calculation of the kinetic viscosity, we obtain 
the VxVy correlation factor 



c(7Ve) = b(i - + ^ + 0{N-^). (27) x^,cc. x,, cc, and % 



FIG. 3: (Color online) Dependence of the viscosity in MPC- 
AT±a on At* in two- or three-dimensional space for (a) n = 5 
and (b) n = 1. Symbols represent the numerical data of 
MPC-AT-ha in 2D (o. A) or 3D (•, T) and MPC-AT-a in 
2D (□, o) or 3D (x, respectively. Solid and dashed lines 
represent analytical results for MPC-AT-|-a and MPC-AT-a, 
respectively. Error bars are smaller than the size of symbols. 



where the numerator and denominator are averaged over 



The kinetic viscosity rykin for large n is then given by 
Eqs. dni) and ^ with Cm = c(n). For MFC- AT -fa and 
MPC-LD+a, this implies for large n that 



independently. When {uS) is also 
pre-averaged over yi^cc, (w) = 7/2 is obtained. How- 
ever, Eq. ((2T|) together with Eq. (|30p contains an integral 
with Hi cc^, which yields an additional correction term of 

0(A^e-'), 



AT+a 
Vkin 



nksTAt r n l^ 

1/ In- (d + 2)/4 ~ 2/' 
^2 , 



(28) 



'?k: 



LD+a _ nkBT r mn(l + 7Ai/2m)77 At^ 
~ 1/ l2n-d- l + d7Ai/4m ^ ~2"/ 



dy^,cc = + ^) + 0{N,~'). (32) 



(29) 



Note that rj and D of MPC-LD±a have a different depen- 
dence on the time step At than other MPC algorithms, 
since their correlation factors A and B depend on At (see 
Table HI). 

The mean velocity difference for MPC+a is given by 



Then, the coUisional viscosity t^coI of MPC+a for large n 
is given by 



?7col 



Am{n — 7/5) 
2Al/-^At 



(33) 



{vf7 



-(l-^M(7 



(30) 



The z component of the velocity is pre-averaged, the an- 
gular velocity is in the vorticity direction, w = we^, and 

(vA = ^y. rre-r.. SO that 



^Uj, CC^X ; 



3 -j.c 
2 



(A^e - i);cVi2 . 



+ (2iV, - l)/e712 



(31) 



Next, we derive the correction terms for small n. For 
A^c = 1 or 2, Eqs. ^ and ^ do not give the cor- 
rect correlation factors s(A'c) and c{Nc) for MPC+a — 
unlike for MPC— a. First, there is no velocity trans- 
fer for A^c = 1, i-e. s(l) = c(l) = 0. Second, in 
the energy-conserving versions of MPC (MPC-SR+a and 
MPC-RA+a), all dNc degrees of freedom are determined 
for A^c = 2 by the conservation of energy (one degree 
of freedom), and translational {d degrees) and angu- 
lar (d — 1 degrees) momentum, so that s(2) = c(2) = 
0. In the MPC versions with an intrinsic thermostat 
(MPC- AT+a and MPC-LD+a), one degree of freedom 
remains for the velocity transfer for A^c = 2, so that 
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FIG. 4: (Color online) Dependence in MPC-AT±a of the vis- 
cosities on the particle number per cell, n. (a),(c) r/kin and 
(b) 7)coi- Symbols represent the numerical data of MPC-AT+a 
(o, A) and MPC-AT-a (□, o) for At* = 0.1 and At* = 1 
in 2D, respectively, and the numerical data of MPC-AT+a 
at At* = 1 (x) in 3D. In (a),(b), the viscosity is rescaled 
by At* and 1/At*, respectively, in order to facilitate a pre- 
sentation of data for different At* on the same scale. Solid 
and dashed lines in (a) and (b) represent analytical results for 
MPC-AT-l-a and MPC-AT-a, respectively. Solid and dashed 
lines in (c) represent analytical results with or without the 
correction term hm, respectively. Error bars are smaller than 
the size of symbols. 




FIG. 5; (Color online) Viscosity of MPC-LD±a as a function 
of (a) 7* and (b), (c) n in three-dimensional space at (a) 
n = 3 and At* = 0.1 and (b),(c) 7* = 1. Symbols represent 
the numerical data of MPC-LD-|-a (o. A) and MPC-LD-a 
(□, o). Dashed and solid lines represent analytical results 
for (a) r/kin or r^coi, and (b), (c) At* = 0.1 or 1, respectively. 
Error bars are smaller than the size of symbols. 



Am, 



t „/7 2n 3n^\ 



?7col 



s(2) A/2d and c(2) ^ {A + B/d)/{d + 2). Thus, 
Sin = Sfc^3 ^'^^ energy-conserving versions 
of MFC, and s,n = p[2)A/dn + Y.V=^ s(k)P{k)k/n for 
MFC versions with an intrinsic thermostat. For MFC- 
SR-|-a and MFC-RA-|-a, the diffusion constant D, and 
the viscosities 7715111 and 77C0I are given by Eqs. (|23p and 
dni) with 



For MFC-AT+a and MFC-LD+a, the diffusion constant 
D and the viscosity contributions T^kin and t^coI are given 
by Eqs. ^ and ^ with 




?/col 



Ail 



1 



2n 
{d-l){d 



(35) 



2)n 



2d 



B{l-e-^\l + n)) + (A+^y 

^{Ad-^^^} 
Am 




d + 2 

B{3d + 2) 1 - e-"(l + n + t? 12) 



(36) 



24:l/-'^At 



(37) 
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FIG. 6: (Color online) Viscosity ri in two-dimensional space 
as a function of At* for MPC-SR±a with 6 = tv/2, MPC- 
RA±a with 9o = tt, and MPC-DR with n = 1 or n = 5. 
Solid and dashed lines represent analytical results for n — 5 
and n — 1, respectively, (a) Symbols represent the numerical 
data of MPC-RA+a (A, A,+), MPC-DR (•,o,<), MPC-RA-a 
(o,x), MPC-SR-a (□), and MPC-SR+a (■). Error bars are 
smaller than the size of symbols. 



l5 



2n 
T 




For MPC-DR, the rotation angle 9 is uniformly dis- 
tributed in — TT < 6* < TT under the molecular-chaos as- 
sumption. Thus, the transport coefficients of MPC-DR 
coincide with those of MPC-RA+a at 6*0 = tt. Thus, 
the diffusion constant D, and the viscosities 77kin and Tycoi 
of MPC-DR are given by Eqs. and ([Ml) with 

A = B = I. Here, the term c,„ can be written in a 
simpler form, = {n — 1 + e^"(l — n^/2)}/n. 



C. Numerical Results 

Figs. [SHS] show the viscosities rykin and 77C0I for five 
MPC methods with or without the angular-momentum 
conservation. The results are displayed in form of di- 




FIG. 7; (Color online) Viscosity difference A?7coi ~ Vcoi — Vcoi 
of MPC-AT— a in two and three dimensions. Symbols with 
dashed or solid lines represent the numerical data in 2D (o, 
□ ) and 3D (o, A) at n = 1 or n = 5, respectively. The inset 
shows the dependence of A?7coi on the average particle number 
n per cell, for At* =0.1 (□) and At* = 1 (o). 
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<1 
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FIG. 8: (Color online) Dependence of the diffusion constant 
D of MPC-AT±a on At at n = 5 in three dimensions. Sym- 
bols and lines represent numerical and analytical data, re- 
spectively. Error bars are smaller than the size of symbols. 



mensionless quantities with length and time units Ic and 
To ~ lc\/ m/k^T, respectively. The main parameters 
which control the properties of MPC fluids, the time 
step and friction constant, have the dimensionless form 
At* — At/ro and 7* = ^Tf)/m. Similarly, the viscosity 
and diffusion constant of a particle are shown in units of 
?7o = \fmk^ jl^^^ and Dq = lc^Jk^T/■m, respectively. 
The error bars of the simulation results are estimated 
from three independent runs. 

Analytical results are calculated from Eqs. ((23l) and 
(fTT]) together with Eq. or from Eqs. ([351) to jSTl), and 
show generally good agreement with the numerical data, 
in particular for At ~ 1 and large n. For smaller time 
step At* =0.1, the most significant deviations between 
numerical and analytical results are found for the kinetic 
viscosity ?7kin, both for MPC- AT-a and MPC-AT-|-a, as 
shown in Fig. Ufa). Similar deviations between analyti- 
cal and numerical results for Tytin have been observed for 
DPD in Refs. [13, [lH, and have been explained by corre- 
lation effects between collisions [13]. At At* = 0.1, a pair 
of particles can collide sequentially several times; in par- 
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ticular for n < 1, pairwise collision occur frequently with- 
out involving any other particles. Thus, the molecular- 
chaos assumption is weakly violated. There are also de- 
viations between analytical and numerical results for the 
viscosity difference ^coi — Vcoi of MPC— a at small At or 
small n (see Fig. [T])- This is also caused by a violation of 
the molecular-chaos assumption. 



Angular-momentum conservation does not affect the 
kinetic viscosity 7/kin of MPC- AT in 2D at large n, com- 
pare Eqs. (|19p and (|28)) . Numerical results are shown 
in Figs. EKa) and |D(a). However, the correction term 
in Eq. ([55]) predicts a small difference of jykin for MPC- 
AT-a and MPC-AT+a for small n ~ 1, see Figs. M.'b) 
and [IK a). The viscosities Tykin and rjcoi of MPC- AT— a, 
and 77coi of MPC-AT+a for large n, show no dependence 
on the space dimension d (except for the scale factor 
therefore, the corresponding symbols and lines in Fig. [3] 
coincide. 



In two dimensions, MPC-SR with 9 and MPC- 

RA with 9o — are characterized hy A — 1, and by 
B = 2 and B = 1, respectively. Thus, they have the same 
collisional viscosity rjcoi for both their '—a' and '-|-a' ver- 
sion, but a different kinetic viscosity Tykim see Fig. [HI Al- 
though MPC-DR has the same viscosity of MPC-RA-|-a 
theoretically, the numerical data of MPC-DR shown in 
Fig. [3] display a slightly larger deviation from the the- 
oretical results for Tycoi and a smaller deviation for ?7kin 
than the data of MPC-RA+a. 



Eq. (jl7p together with (|36p predicts a minimum of T^kin 
around rt = 1, as shown in Figs, djc) and [5Kb). How- 
ever, this minimum is not seen in numerical data and 
could be caused by the negligence of higher-order terms 
in Eq. ([77)1. We therefore investigate the dependence of 
the next-order term h/Nc^^ where /i is a free parameter. 
The average is estimated by /im = X]fc°=3 ^(^)^/^^ — 
{1 - e""(l +n + n'^ + ri^/G - n'^/72)}h/n^, which yields 
the asymptotic dependence hm — hn^ / 18 for small n and 
= h/n'^ for n —> oo. The correction term is then 
added to Eq. ([55]) with /i as a fit parameter. Fig. [5Kc) 
shows that this correction term with h = —0.6 in 2D and 
/i = — 1 in 3D removes the minimum and gives better 
agreement with the numerical data of MPC-AT+a. 



Fig. [8| shows the self-diffusion constant D of MPC- 
AT+a. The '+a' fluid displays faster diffusion than the 
'—a' fluid. The diffusion constant D is numerically cal- 
culated from the mean square displacement of a particle, 
({ri(i) — ri(0)}^) = 2dDt, in a cubic simulation box with 
side length L = 201^- Deviations from the analytical re- 
sults calculated with the molecular-chaos assumption are 
seen for small At* . 



III. DISSIPATIVE PARTICLE DYNAMICS 
(DPD) 

A. Simulation Method 

The DPD thermostat is a modified Langevin thermo- 
stat, where friction and noise forces are applied to the rel- 
ative velocities of pairs of neighboring particles d, H, 0] ■ 
The equation of motion for the i-th. particle with mass m 
is given by 

dVi dU sr-^ ( „ , / Nl - 

(38) 

where = - v^, = r,; - r^, and f,y = 
Tij/rij, with weight Wij = 'w{rij). The Gaussian white 
noise S,ij{t) obeys the fluctuation-dissipation theorem, 
with {^,j{t)) = and {i,j{t%,f{t')) = 2kBT{5u'6ry + 
5iji5iji)6{t — t'). This thermostat is applied only in the 
direction fy to conserve the angular momentum. We 
denote this original method here DPD+a. 

In DPD, a linear weight function = Wi{rij) = 
7(1 — Ty /rcut) is typically employed, which vanishes be- 
yond the cutoff distance — Vcut- Furthermore, DPD is 
usually combined with a soft repulsive potential U; how- 
ever, we only consider the ideal-gas equation state (with 
potential U — 0) in this paper. 

The DPD equation ([55)1 is discretized by the Shard- 
low's SI splitting algorithm Q, where each thermostat 
of the ij pair is integrated separately, 

v"™ = V, + {-adp(ry )vij • + 6dp(»'ii)Cij,n}fii, 

(39) 

with 

_Wj^Atlm_ ^ _ y/w.,jAt/m 

fldp {nj }- ^^^^^^^1^^ OdP j - ^ ^ ^^^^ . 

(40) 

The discretized Gaussian noise ^ij^n is determined by 
the variance {^ij^n^i'j',n'} — 2k^T{^6ii'djj' ~\~ Sij'Sij')6nn' • 
This splitting algorithm belongs to the class of gen- 
eralized Lowe- Anderson thermostats [l^, because the 
factors ttdpirij) and bdp{rij) satisfy the relation 6dp = 
^adp(l - adp)/m 

DPD can be modified to remove angular-momentum 
conservation. We denoted this technique here DPD— a. 
It has been introduced in Ref. [1^ to explore the similar- 
ities and differences between DPD and MPC methods. 
In this case, the equation of motion reads 

dvi dU \ - r . , , , 

The splitting algorithm can also be applied to DPD— a 
as vf°^ = V, - adp(ry)v,j + 6dp(r,j )!„■„. 
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The combination of DPD+a and DPD— a, denoted 
'transverse DPD', with an equation of motion determined 
by the difference of the right-hand sides of Eqs. ((4T|l 
and has been suggested very recently [i^. A sim- 
ilar anisotropic friction has been used in the standard 
Langevin equation to tre at p olymer entanglement im- 
plicitly in polymer melts [50| and dilute polymer solu- 
tions [51]. 

The DPD thermostat can be generalized into a multi- 
body thermostat (denoted DPD-MT-a) which is 
defined by the equation of motion 



I—— 

dt 



dU_ 



+ (42) 



where un 



and vp = ^j^iWij\j/w^ is the 
weighted mean velocity. The second term on the right- 
hand side of Eq. ([55]) is the friction term between the 
j-th particle and its neighbors, and thermostats in 

Eq. ([51]) are unified into a single thermostat, where A^„b is 
the average number of the neighbors with rij < rent- The 



third and fourth terms on the right-hand side of Eq. 
are needed to conserve the translational momentum. 

Angular momentum can be conserved in DPD-MT, 
when the thermostat for the z-th particle is applied only 
in the direction r^^c = — rp, where the weighted center 
of mass is rp = '^j^i '^ij'^j/'^^- The equation of motion 
of DPD-MT-Fa is thus given by 



dU 



(43) 



Shardlow's SI splitting algorithm Q can be apphed to 
both DPD-MT-a and DPD-MT-ha. Eq. gj) of DPD- 
MT— a is discretized such that each thermostat of the 
1,1*^ pair is integrated separately, 

v-* = V, - af (v, - vP) + (44) 
= v, + ^{«f(v.-vP)-5fe,„}- 
The factors a™* and 6™' are given by 
„mt _ w°At/m 



, (45) 



1 + z/iu;0Ai/2m' ~ 1 -f i^,w^At/2m 
where 1/^ = 1 + J^j^i "^ii^ I ^^Xf ■ 

B. Transport Coefficients 



We now derive analytical expressions for the viscosity 
T] and self-diffusion constant D of DPD— a and DPD- 
MT±a with ideal-gas equation of state (with potential 



[/ = 0). The corresponding derivations for DPD-|-a (2l| 
can be straightforwardly carried over to this case. 

The correlations of DPD±a results from a multitude 
of pairwise collisions, so that 1 — s 

Cm — (n^C' 



— (IljSjj) and 1 



jL-ij/. Eq. ((39)) together with a molecular-chaos 



assumption implies Sij = 1 — adpif 



= 1 - adp(x| 



ai^) + Aadp^xfjijfj for DPD+a, and Sij = 1 - a^p, Cij = 
1 — 2adp-|-2adp^ for DPD— a. For an ideal gas, the number 
of particles k per volume AV is given by the Poisson 
distribution, P{k) = e-"^'^(nAF)'=/fc!, so that (c'^) = 
exp{(— 1 -|- c)r7,AV"} for some constant c. This implies 

1 - S„i = CXP(-1 +^^.(s,y)). 

The collisional stress cr™/ is the momentum flux due to 
DPD collisions crossing a plane at y = yo- After inter- 
change of the order of integration, tr™' is given by 



^col 



2At 



(46) 



where Eq. ()39p and (vij^x) = iVij have been used. Thus, 
the diffusion constant and viscosity of DPD-|-a are given 
by Eq. ^ with [U 



D 



Vcol 



ksTAt 



1 — exp ■ 



2d{d + 2) 



exp(-n[adp(r)]g/d) 
2adp(r) 4adp(r)2 



wr 



did + 2) 



1 -|- wAt/m 
g(r)w{r)dV. 



,(47) 
(48) 

(49) 
(50) 



Similarly, for DPD— a, the viscosity and diffusion con- 
stant are found to be 



D 



ksTAt 
m 

1 — exp 



1 



1 — exp( 



n[adpir)]g) 
|2n [-adp(r) +adp(0^]g} 



??col = 773 



n 
Yd 



1 -|- wAt/m 



(51) 
(52) 
(53) 



The only differences between the expressions for D, Cm 
and r^coi in DPD-|-a and DPD— a are prefactors containing 
d and d -t- 2. 

To simplify the equations of DPD-MT, the factors a™* 
and i^i are pre-averaged as 

r,„2i 



n[w]gAt/m 
1 + i^,nn[w]gAt/2r 



1 



Then D and rjcoi of DPD-MT+a are given by 



D 



ksTAt 



1 — exp 



1 - exp{-u„,anjd) 
2a,nZ^m . 2a 



m 



d{d + 2) 



Vcol 



d{d + 2)[w]g{l + iy^n[w]gAt/2m) 



(54) 

(55) 
(56) 
(57) 
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FIG. 9: (Color online) Dependence of the viscosity 77 on At* 
of DPD±a and DPD-MT±a in three-dimensional space for 
nVcut = 3 and jTo/m = 9. Symbols represent the numerical 
data of DPD+a (o), DPD-a (x, □), DPD-MT+a (A), and 
DPD— a (o). Dashed and solid lines represent analytical re- 
sults for DPD— a and other DPD methods, respectively. Error 
bars are smaller than the size of symbols. 



Finally, for DPD-MT-a, we find 
kBTAt ( 1 



D 



Vcol 



m il- exp(-zy,„am) 

1 - exp(-2ai„t'm + flin^t'm^), 

n[w'^r'^]g 
d[w]g{l + u^n[w]gAt/2m) 



C. Numerical Results 



(58) 
(59) 
(60) 



Fig. [9] shows the viscosity of various DPD fluids witli 
an ideal-gas equation of state and the linear weight 
wi{rij). The viscosity and time step are normalized by 
770 = y/ mkBT / rcut'^~^ and tq = Tcut^/mjk^ , respec- 
tively. The dimensionless time step is At* = Ai/ro, as 
before. There is in general good agreement between ana- 
lytical and numerical results. However, small deviations 
are visible. One reason for these deviations is that the 
molecular-chaos assumption is not perfectly valid p^ . 
In the case of DPD-MT±a, another reason is the pre- 
averaging procedure used in the derivation of the analyt- 
ical expressions, which neglects some correlations. 

The kinetic (collisional) viscosities of DPD-|-a and 
DPD-MT-|-a are larger (smaller) than those of the '—a' 
versions, since angular-momentum conservation reduces 



the momentum transfer in DPD collisions. A similar be- 
havior has also been found for MPC±a in Sec. HT] 



IV. THERMOSTATING MESOSCALE FLUIDS 
UNDER FLOW 



In experiments, systems are usually thermostated on 
their boundaries. However, in simulations, thermostats 
typically act on all fluid particles in order to avoid tem- 
perature gradients. In flows, the temperature is defined 
under the assumption of local equilibrium. In the MPC 
and DPD families, the length scales which define this "lo- 
cal" environment are Ic and Tcut, respectively. On these 
scales, the thermal fluctuations should be separated from 
the macroscopic flow, and the thermostats should act on 
the local kinetic energy to fix the temperature. 

The conditions on the shear rate 7 for this local equi- 
librium to hold are obtained as follows. All of ther- 
mostats of the MPC family are profile-unbiased ther- 
mostats (PUT) [131 . Thus, the condition for a max- 
imum shear rate of PUT [531 also apply to MPC. In 
simple shear flow with low Reynolds number, the par- 
ticle velocities are characterized by (vi(ri)) — jyiex and 
(v.,(ri)2) = dkBT/m + j'^yj^. In MPC, the particle veloc- 
ity Vi^c relative to the center-of-mass velocity of a MPC 
collision cell is employed to calculate the kinetic energy 
in the local rest frame. 



1 



j^^ccll 



m 
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(61) 



where the average is taken over all particles in a cell. 
For ^ ^k^T /m, the second term in Eq. (|6ip is neg- 
ligible, and the thermal fluctuations and shear are well 
separated. On the other hand, for -^l^ > y/k^r/rn, the 
thermostats couple with the macroscopic flow and may 
modify the flow behavior. 

In DPD-f-a, the relative velocity v,j of neighboring par- 
ticles is employed instead, 



((V,J •fy)2) 



m 



\w\ 



(62) 



For the linear weight wi{rij) and uniform radial dis- 



tribution function g{r), the second term in Eq. (|62[) is 
{{d + l)/{d + 2f{d + 3)}(7rcut)^- Thus, the condition 
for thermostats to provide local equilibrium conditions is 
7'"cut < y/k^TJm. 

To study the hydrodynamic behavior of complex fluids, 
the parameter ranges of simulations should of course also 
match physical conditions of experiments. Thus, the sim- 
ulation parameters have to be chosen such as to adjust di- 
mensionless hydrodynamic quantities, like the Reynolds 
number, the Schmidt number, and the Knudsen number. 
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V. SUMMARY 

MPC and DPD are very versatile simulation tech- 
niques for mcsoscalc hydrodynamics. By employing dif- 
ferent types of collision rules and thermostats, it is pos- 
sible to construct a variety of algorithms with different 
properies. One of the important properties is whether an 
algorithm does or does not conserve angular momentum. 
The angular momentum conservation can be switched on 
or off in each variant of MPC and DPD. 

In addition to MPC algorithms suggested previously, 
we have introduced here an angular-momentum con- 
serving version of the widely used stochastic-rotation- 
dynamics algorithm of MPC. This algorithm has to be 
used with some caution, because compared to other 
MPC-|-a techniques, it does not give a uniform radial 
distribution function. However, the deviations are small 
for sufficiently large particle numbers per cell and not too 



small time step. 

We have derived analytical expressions for the viscos- 
ity 77 and the self-diffusion constant D of various MPC 
and DPD methods. The theoretical results show very 
good agreement with numerical results. Many similari- 
ties between MPC and DPD are seen in the derivation of 
T) and D and the relation between the '—a' and '+a' ver- 
sions. We believe that these similarities apply generally 
for particle-based hydrodynamics methods. 
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